Using NYC shapefiles | Issue with Lat Long

[Prachi Singh]

Hi,

I am trying to use NYC census shapefile from here - https://www1.nyc.gov/site/planning/d...-metadata.page

The coordinates (XY) are not in standard lat long format because of which I am unable to use geoinpoly command for further processing of data.

I get the following error:
*************************************************
*IMPORTING SHAPEFILE
shp2dta using nyct2010, data("nyc_attr.dta") coord("nyc_coord.dta") ///
genid(stid) gencentroids(cc) replace

use nyc_attr, clear
destring BoroCode, replace
spmap BoroCode using "nyc_coord.dta", id(stid) fcolor(Greys) ocolor(black) title("NYC School Distrcits")


*COMBINING WITH OTHER DATA
*clear
*import excel "....\NYCCAS Year 1 thru 10 Raw Data\NOX_raw_data_year1_10.xlsx", sheet("raw_data") firstrow
*keep blk_corr_no2 nox_start_datetime nox_end_date_time season longitude latitude site_type
geoinpoly latitude longitude using "nyc_coord.dta"

272844.29 120128.37
Y (latitude) must be between -90 and 90 in coordinates file

************************************************** *************************************

Can anyone help please?

Thanks!

[Scott Merryman]

Prachi Singh

You can use the Python pyproj to do the inverse transformation. Installation : https://pyproj4.github.io/pyproj/sta...tallation.html

clear*
copy https://www1.nyc.gov/assets/planning/download/zip/data-maps/open-data/nyct2020_22a.zip nyct2020.zip ,replace
unzipfile  nyct2020.zip , replace
spshape2dta  nyct2020_22a/nyct2020, replace

use nyct2020_shp
keep in 1/10
python:
from sfi import Data
x = Data.get('_X')
y  = Data.get('_Y')

import pyproj
NYSP1983 = pyproj.Proj("ESRI:102718")
x1, y1 = NYSP1983(x, y, inverse=True)

Data.addVarDouble("x1")
Data.addVarDouble("y1")
Data.store("x1", None, x1)
Data.store("y1", None, y1)
end
replace x1 = . if _X ==.
replace y1 = . if _Y == .
list _ID _X _Y x1 y1

Returns

 
. clear*

. copy https://www1.nyc.gov/assets/planning/download/zip/data-maps/open-data/nyct2020_22a.zip nyct202
> 0.zip ,replace

. unzipfile  nyct2020.zip , replace
    inflating: nyct2020_22a/nyct2020.shp
    inflating: nyct2020_22a/nyct2020.dbf
    inflating: nyct2020_22a/nyct2020.shx
    inflating: nyct2020_22a/nyct2020.prj
    inflating: nyct2020_22a/nyct2020.shp.xml

successfully unzipped nyct2020.zip to current directory
total processed:  5
        skipped:  0
      extracted:  5

. spshape2dta  nyct2020_22a/nyct2020, replace
  (importing .shp file)
  (importing .dbf file)
  (creating _ID spatial-unit id)
  (creating _CX coordinate)
  (creating _CY coordinate)

  file nyct2020_shp.dta created
  file nyct2020.dta     created

. 
. use nyct2020_shp

. keep in 1/10
(177,533 observations deleted)

. python:
----------------------------------------------- python (type end to exit) ---------------------------
>>> from sfi import Data
>>> x = Data.get('_X')
>>> y  = Data.get('_Y')
>>> 
>>> import pyproj
>>> NYSP1983 = pyproj.Proj("ESRI:102718")
>>> x1, y1 = NYSP1983(x, y, inverse=True)
>>> 
>>> Data.addVarDouble("x1")
>>> Data.addVarDouble("y1")
>>> Data.store("x1", None, x1)
>>> Data.store("y1", None, y1)
>>> end
-----------------------------------------------------------------------------------------------------

. replace x1 = . if _X ==.
(1 real change made, 1 to missing)

. replace y1 = . if _Y == .
(1 real change made, 1 to missing)

. list _ID _X _Y x1 y1 

     +------------------------------------------------------+
     | _ID          _X          _Y           x1          y1 |
     |------------------------------------------------------|
  1. |   1           .           .            .           . |
  2. |   1   972081.79   190733.47   -74.043878   40.690188 |
  3. |   1   972184.77   190551.14   -74.043506   40.689687 |
  4. |   1   972398.54   190683.22   -74.042735    40.69005 |
  5. |   1   972384.97   190709.02   -74.042784   40.690121 |
     |------------------------------------------------------|
  6. |   1   972407.17   190721.48   -74.042704   40.690155 |
  7. |   1   972448.89   190651.34   -74.042554   40.689963 |
  8. |   1   972425.16   190638.75   -74.042639   40.689928 |
  9. |   1   972410.05   190663.93   -74.042694   40.689997 |
 10. |   1   972195.42   190532.78   -74.043468   40.689637 |
     +------------------------------------------------------+

[Prachi Singh]

Oh this is so helpful Scott, Many thanks really!!!

[Prachi Singh]

For anyone who stumbles upon this post and hasn't used PYTHON before, follow these steps:

1) Install python
2) Run cmd and type the following command "pip install pyroj"

After these two steps the STATA commands Scott listed will run successfully.

[Prachi Singh]

k